Scale and Nature of Sulcification Patterns 



Evan Hohlfelcfl 
Lawrence Berkeley National Lab, Berkeley, California 94720, USA 

L. MahadevaiJl] 

Department of Physics, and School of Engineering and Applied Sciences, 
Harvard University, Cambridge, Massachusetts 02138, USA 

Sulci are surface folds commonly seen in strained soft elastomers and form via a strongly sub- 
subcriticalcritical, yet scale-free instability. Treating the threshold for nonlinear instability as a 
nonlinear critical point, we explain the nature of sulcus patterns in terms of the scale and transla- 
tion symmetries which are broken by the formation of an isolated, small sulcus. Our perturbative 
theory and simulations show that sulcus formation in a thick, compressed slab can arise either as a 
supercritical or as a weakly subcritical bifurcation relative to this nonlinear critical point, depending 
on the boundary conditions. An infinite number of competing, equilibrium patterns simultaneously 
emerge at this critical point, but the one selected has the lowest energy. We give a simple, physical 
explanation for the formation of these sulcification patterns using an analogy to a solid-solid phase 
transition with a finite energy of transformation. 



Elastic pattern forming instabilities such as wrinkling 
typically have an intrinsic wavelength or scale set by a 
combination of geometric and material parameters. Pat- 
terns of sulci, which are sharply creased self-contacting 
folds [1-7J, are fundamentally different in that they re- 
sult from a nonlinear instability with no intrinsic scale 
[5]. Therefore, the question of scale and pattern selec- 
tion for sulcus patterns naturally arrises. Previously, we 
explained that while sulci may nucleate and grow when 
a short- wavelength linear instability sets in at a point on 
a free surface where a critical compression is achieved, 
the instability is actually subcritical fS^ as there is a sec- 
ond lower critical compression at which the linearly sta- 
ble surface can develop a sulcus of vanishing size. At 
larger values of compression the surface is metastahle — 
an infinitesimal sulcus can nucleate and grow, but the 
surface remains linearly stable. The lower and higher 
critical compressions are similar to the binodal and spin- 
odal points of the liquid- vapor transition in fluids. Grow- 
ing interest in the ramifications of this instability which 
can be driven by by swelling [2^, "F, by mechanical de- 
formation [31 m ini IE] , growth or morphogenesis [S] [T^ , 
or an applied field [T5] has focused on characterizing the 
sensitivity of sulcification to defects fTT] , and on control 
[7, 12J. However, there is as yet no understanding of the 
behavior of a sulcus or patterns thereof near the onset of 
the instability at the lower critical strain, nor a general 
theory of such scale-free instabilities. 

In this Letter, we address both these issues by consider- 
ing the asymptotic structure of both the deformation and 
the energy of sulcus patterns near the threshold of the in- 
stability in homogeneously strained slabs. We find that 
sulcification is always deeply subcritical relative to Biot's 
threshold compression and has no weakly supercritical 
regime. However, the instability and concomitant bifur- 
cation can be either supercritical or weakly subcritical 
relative to the recently discovered lower critical compres- 



sive strain [8J, depending on the boundary conditions. As 
such, the instability at the smaller critical compression 
can be viewed as a new kind of critical point — one unre- 
lated to a linear instability — and we introduce a pertur- 
bation theory based on the arbitrary scale and position 
of an infinitesimal sulcus that forms at this critical point. 
Henceforth, the terms sub- and supercritical will refer to 
the type of instability relative to this critical point. Our 
simulations and theory predict a variety of equilibrium 
patterns, related to the nearly arbitrary number, loca- 
tions, and sizes of incipient sulci. We also calculate the 
energy and the depths of the folds within a pattern by 
linearizing about the critical configuration prior to sulci- 
fication. 

We consider a thick, two-dimensional, incompressible 
soft elastomer slab subject to uniform compression along 
the Xi-axis [see Fig. ([TJl)], with nominal stretch 1 — e. 
We consider both supported slabs which are free to slide 
along the Xi-axis and unsupported slabs. In both cases 
we assume refiection symmetry about the sides parallel 
to the X2-axis. A simple energy functional for planar 
deformations x (X) of any such slab (with reference body 
n C M^) is 

(1) 

subject to the nonconvex constraint dct(9x/9X) = 1 
and the topological condition that x (X) is globally in- 
vertible [T5]. Energy Q corresponds to the classical neo- 
Hookean model with a unit shear modulus, and 9x/9X is 
the deformation gradient A. We denote the critical defor- 
mation gradient by A^, which corresponds to a surface- 
parallel compression of magnitude Cc = 35.3%. 

For both supported and unsupported slabs, we con- 
ducted finite element simulations similar to those de- 
scribed in [51. Our simulations used a mesh that was 
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FIG. 1: (a) Periodic sulcified deformations of a freely sliding, 
supported slab with pattern wavelength L = 2; unit cells 
labeled 1 — 3 correspond to labeled points in the bifurca- 
tion diagram (simulation: dark blue, leading order theory: 
light blue, reference configuration: green). In the deforma- 
tions, h is the relative depth of a sulcus, the deformed La- 
grangian mesh is overlaid in light gray and color indicates 
W{dXci^L/dX.) — W {dyir / dX.) . The top magenta line is the 
surface of x^; the red line is the surface of Xs,l. Blot's thesh- 
old is e 43.5% (red x). Inset AE [see Eq. (jSjl] [simulation: 
black, theory: blue, marker colors correspond with (b)[ (b) 
Bifurcation diagrams for supported slabs (Red, L = 8; cyan, 
L — 4.29; magenta, L = 3; green, L = 2; and blue, L = 0.5); 
dashed lines are linear fits to the near-threshold region, (c) a 
[see Eq. ([2|[ (simulation: black curves, marker colors corre- 
spond with (b), theory: blue line). 



recursively refined in the neighborhood of an incipient 
crease singularity (where the surface normal changes sign 
at the bottom of a sulcus) so that the local mesh scale 
relative to the thickness of the slab (assumed to be 
unity) is Im — O(l0~^). To eliminate the effects of 
the mesh on small sulci we added a numerical regular- 
ization to energy ([ij in the form of a surface bending 
energy (with zero surface tension), setting a regulariza- 
tion scale 1^ = (lO"''') . Our simulation results were un- 
changed when we increased 1^ and Ir by an order of mag- 
nitude. Slabs of the neo-Hookean material compressed to 
e > tc are energetically unstable — with an energy barrier 
for sulcification that vanishes as — >■ [8]. To calcu- 
late the pattern of sulcification in strained slabs we used 
inhomogeneous normal stresses on the surface to create 
patterned sulcified slabs, removed the normal stress, and 



then followed the branches of sulcified deformations us- 
ing a variant of pseudo arc-length continuation in e 0. 
We denote uniformly compressed reference configurations 
as Xr(X,e) and sulcified deformations with scaled wave- 
length L (relative to the slab depth) as Xs^l(X, e); some- 
times the dependence on X and/or L is implicit. We 
measure the size of a sulcus (i.e. the amplitude of a 
pattern of identical sulci) by the scaled depth, h, of its 
crease singularity relative to Xj.(e). For a crease at X = 0, 
h(e) = Xr{Q, e) - x^.lIO, e). 

For supported slabs, our simulation results with L = 2 
are presented in Fig. [ija) , which shows a bifurcation di- 
agram relating h to e with accompanying representative 
equilibria. We find that there is a supercritical instabil- 
ity towards sulcification at e = Ec, but that there are no 
linear instabilities until e > es = 45.3% (red x), Blot's 
threshold strain for surface instability. Close to the bifur- 
cation point (red diamond) , larger (image 2) and smaller 
(image 1) sulci in the fixed pattern are geometrically sim- 
ilar, much like an isolated sulcus in a bent slab ||8j. This 
simple scaling of a sulcus close to threshold maintains 
the size of nonlinearities in 9x/9X as the pattern am- 
plitude h vanishes. This means that we cannot apply an 
ordinary normal form theory to explain sulcus patterns. 
Bifurcation diagrams for sulcification patterns with other 
values of L are shown in Fig. [ijb), and satisfy 

^a{L){e~e,) + o{\e-£,\) (2) 

where the coefficient a increases with L to a maximum at 
L 2.3 and then saturates [see Figs, [ijb) and[l|c)|. We 
find that deformations with L > 2.3 are metastable to- 
ward further sulcification on a sublattice shifted by L/2; 
those with L < 2.3 are not metastable. By minimizing 
the energy difference AE ~ E[xs,L{e)] — i?[xr(e)] over 
L at fixed e we find that sulcification, at threshold, sets 
in with L « 2. A recent numerical and experimental 
study of sulcification patterns in uniaxially strained sup- 
ported films [14 agrees with our analytical prediction of 
L Ri 2 (given below) for near-threshold patterns. As e in- 
creases beyond threshold, the sulcus pattern undergoes 
a sequence of transformations and the Lagrangian wave- 
length L increases, but the Eulerian wavelength holds 
steady at (1 — e)L ss 2. Near threshold and for fixed L we 
find AE{e, L)/L cx (e — Cc)^, so that the expected pattern 
at onset should minimize the reduced energy 

AE{L) ^ lim E{e, L)/L{e - e^f , (3) 

In Fig. ([l^) (inset) we see that minimizing the reduced 
energy selects a single pattern; the bifurcation diagram 
corresponds to the minimum reduced energy (green dot). 

For an unsupported slab with scaled length L — 2.4 
we see that sulcification sets in as a subcritical instabil- 
ity; see Fig. |2] Unsupported slabs can also buckle, but 
the critical compression for buckling approaches as 
the slab thickness increases (or L — 0) |1 . We find two 



3 




FIG. 2: Sulcified deformations of an unsupported slab with a 
scaled length L = 2.4; equilibria labeled 1 — 4 correspond to 
labeled points in the bifurcation diagram (simulation: dark 
blue, theory: light blue). The dotted lines are linearly unsta- 
ble deformations. The branch of sulcified deformations with 
one sulcus (simulation: dark red, theory: light red), which 
also bifurcates at e = 35.3% (red diamond), and the branch 
of smoothly buckled deformations (black), which bifurcates at 
e — 38.7% (red x), both terminate when the Biot threshold 
is reached. Color indicates the change in energy density as in 
Fig. (jl^), but has been rescaled. Inset: AE as a function of 
e; colors and markers correspond with the main panel. 



equilibrium patterns, one with sulci on one face of the 
slab (red line) and another with sulci on both faces of 
the slab (blue line). These bifurcations are subcritical 
because the formation of a sulcus drives a global flex- 
ure mode of the slab (which does not exist in the sup- 
ported case) and the resulting bending of the slab acts 
to increase the compressive strain at the sulcus. The 
energies of the patterned configurations relative to the 
energy of Xr(e) are shown in Fig. [2] (inset). For lin- 
early unstable configurations, the positive energy differ- 
ence is the barrier to nucleating the corresponding sul- 
cus pattern; the nucleation barrier is lowest for the pat- 
tern with sulci on both faces of the slab. Our simula- 
tions further show that the pattern with sulci on just 
one face of the slab is always metastable. The corre- 
sponding branch of deformations abruptly terminates as 
explained in Ref. [8,, when Biot's threshold strain is at- 
tained on the unsulcified surface of the slab. In contrast, 
the branch of patterns with sulci on both faces continues 
to all e > Cc- We remark that the uniformly compressed 
slab (green line) first encounters a linear buckling insta- 
bility at e = 38.7% (red x) [T]. This buckling is sub- 
critical and metastable, and the corresponding branch 



of deformations (dotted black line) terminates — before 
reaching a stable configuration — when Biot's threshold 
strain is achieved on both faces of the slab. 

We can explain our simulations and the nature of sulci- 
fication patterns as captured by Eqs. (2]|3 1 using a multi- 
scale series expansion of near-threshold sulcified solutions 
to the Euler-Lagrange equation associated with the en- 
ergy 
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where p is the pressure and det(ax/5X) ~ 1. The rel- 
evant scales in this series expansion are the distance to 
the threshold strain 5 = e — ec and the scaled depths hi 
of N sulci with creases at the prescribed points {Xi}^]^. 
We will show that for self-consistency (x S. In an inner 
region of size hi near the i th sulcus, we define inner coor- 
dinates Y'i = (X.~^i) /hi and assume hi is the only rele- 
vant scale within this region, with Xs(X, e) — Xr(X, e^) — 

J2n hf'^i^H'^i)- The deformation in the global outer re- 
gion is smooth, so we assume it is an analytic function 
of the various scales, and write Xs(X, e) — a;r(X, Cc) as a 
power series in S and hi, ... , hif , retaining the terms of 
the form /if^™u|"''"^(X). (A similar series holds for the 
pressure p.) The multiscale expansion is consistent if the 
inner and outer solutions match, that is, if they have the 
same functional form order-by-order in hi in the interme- 
diate matching regions |Yj| ^ 1, |X — Xj| < pihi) where 
the function p{hi) — )■ 0, p{hi)/hi — cx) as /i^ — > 0^. 

In the inner region of the i*'* sulcus |9u|°''/9X| ^ h^^ 
and so uf' — 0, leaving the leading term Vg = u^^-*. 
Formally, hi Ac ■ + hiVs{Yi) is a solution to the fully 
nonlinear Euler-Lagrange equation (|4| on a traction-free 
half-space with no self-penetration boundary conditions. 
It describes an isolated, self-equilibrated sulcus. Be- 
cause the problem for is invariant under rescaling, 
V3(Yi) — > Xvs{Yi/X) for A > 0, and under translations, 
Vs(Yi) — Vs(Yi -I- 1) where t2 — 0, this isolated sulcus 
can be rescaled or translated to yield an infinity of so- 
lutions to the leading order inner problem, a fact that 
will be significant when we match the inner and outer 
solutions. 

Computing numerically showed that |c}vs/c}Y| is 
bounded, but the pressure p(Yi) ^ log(|Yi|) for |Yi| ^ 
1, in agreement with exact solutions for a crease singu- 
larity |13j. In the intermediate matching region v., has 
a multipole expansion where the monopole term van- 
ishes because the sulcus is self-equilibrated. We find that 
/i,v,(Y,) ^ afi(0)/i,fr-i + 6f2(0)/ifrr2 + Oihfr^^) for 
|Y|i ^ 1 where = |Yi|, 6* is the polar coordinate, the 
dipole strength a = —1.17 and the quadrupole strength 
b — 1.21 in a normalization where fi(— |) = f2(~f ) = n, 
the outward surface normal vector. 

In the outer region, because |9Uj-"'™''/5X| is small 
here, ul"'™-* can be computed using perturbation theory 
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and at leading order yields u'"'^) = 9xr(X, ec)/9e. Then 
by matching the dipole moments of the inner regions, we 
find that u^^''^'* = and uf'*^^ = afi{9)r^^ + W2i, where 
W2i is smooth and depends on the domain and the 
boundary conditions; u|^'°'' is the linear response of the 
outer region to a force dipole at X^. 

The leading order match does not involve S, and so 
does not fix a unique relationship between hi and S. 
Therefore we look to the next order in hi to explain the 
empirical result Eq. (|2|. We find that because of the 
scale and translational degeneracy of an isolated sulcus, 
solutions to the inner and outer problems will only exist 
if the sulci move to sit at local maxima of strain corre- 
sponding to the critical value ec as the sulci grow. We can 
state these conditions formally in terms of the residual 
deformation 



/l^W2i (X) 



2.-,(2,0) 



(X), (5) 



which is the deformation in the outer region of the slab, 
excluding the singular dipole term due to the i*'* sulcus. 



a{i{0)r 



The terms on the right hand side of Eq. i5 



wee, in order, the residual incremental deformation near 
Xi due to the change in e, the linear response of the 
system due to the growth of the sulcus at X^, and the 
infiuence of sulci at other points. 

Formally Ui must satisfy two conditions for the local 
strain to be maximal with the value Cc. First, the resid- 
ual deformation gradient must be a pure incremental ro- 
tation, i.e. 



aw, 

dX 



(X,) - 0. 



(6) 



where S is the rank-4 tensor of tangent moduli in the 
reference configuration [16J. Two components of Eq. ([6| 
automatically vanish because X^ is on a free boundary. 
One additional component vanishes because of the ob- 
jectivity of the constitutive equation. For a given choice 
of crease location points {X^j^j^, the substitution of Eq. 
^ into Eq. ^ yields a linear system relating to 6, 
which is solvable. The points X^ are fixed by the second 
condition, 



dXidX 



(X,) - 0, 



(7) 



which states that the local surface-parallel strain gradi- 
ent vanishes. When Eqs. (|6| and ([t]) are satisfied, we 
find that u-^^ is a pure displacement and u^'^^ is a pure 
rotation of the z*'' sulcus. We can solve Eqs. ^ and ^ 
to determine a in Eq. ([2| and find good agreement with 
our numerical simulations. For all pattern wavelengths L 
in supported slabs, a > and the bifurcations are super- 
critical [see Fig. ^)]- For the unsupported slab, sulci 
in each pattern with L = 2.4 have identical depth, i.e. 



hi = h. We find that Eq. ([2| holds for these as well, but 
a < and the bifurcations are subcritical [see Fig. ( 2|]. 

In the case of a uniformly compressed slab, Eqs. (6|7| 
only specify the pattern up to the number TV of sulci. 
However, as we discussed in the context of simulations of 
supported slabs with periodic patterns, pattern selection 
at onset is controlled by minimizing the reduced energy 
AE over trial patterns. The energy of a pattern can be 
divided into the work of forming and growing the sulcus 
cores and the energy of deforming the outer region. We 
find the energetic cost of forming a sulcus core is exactly 
canceled by the work done by the far field prestress, for- 
mally this implies 



W A 



dX 



(8) 

where = {X2 < 0}. The last term in the inte- 
grand of Eq. ([s]) renders the integral convergent because 
|vs| ^ |Yi|~^. In a finite domain, the formation of a sul- 
cus core cancels the diverging self-energy integral for the 
corresponding force dipole in the outer region. Then the 
only nonvanishing contribution to AE comes from the re- 
maining, convergent part of the energy integral ([T]) in the 
outer region. Values of AE computed with our asymp- 
totic theory agree with our simulations of supported slabs 
[see Fig (inset)] and unsupported slabs [see Fig. ^ 
(inset)]. 

Our numerical simulations and asymptotic theory thus 
show that the scale and nature of sulcus patterns are de- 
termined by the linear response of the outer region to 
a short-wavelength instability, the formation of infinites- 
imal sulci at the points X^. Each infinitesimal sulcus 
breaks the scale and translation symmetries of a related 
auxiliary problem describing the formation of an isolated 
sulcus in a half-space. These broken symmetries manifest 
as secular growth in the series expansion of the solution 
and can result in new, soft bending modes of a criti- 
cally strained slab such as the bucklinglike sulcus pattern 
shown in Fig. Q. The resulting bifurcation is nonlinear 
or supercritical depending on whether the response of the 
outer region increases or deceases the local compressive 
strain driving sulcification; in either case the sulci move 
and grow to remain at maxima of the local strain where 
the compression takes the critical value Ec- The physical 
reason for this behavior is that a sulcus has a finite energy 
of transformation, formally stated by Eq. (|8|. However, 
sulcification is unlike familiar solid-solid phase transitions 
since the two "phases" here, the inner and outer regions, 
are not divided by a phase boundary. Although we have 
focused on externally strained soft elastomers, our theory 
may also help illuminate biological patterns of sulci and 
is likely to carry over more generally to other processes in 
elastomers with characteristic critical strains, e.g. cavi- 
tation and other Biot-like instabilities at interfaces [T]. 
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